set.seed(1234)
library(raster)
## Warning: package 'raster' was built under R version 4.0.4
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.0.4
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
#library(ggtext)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.5
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ncdf4)
## Warning: package 'ncdf4' was built under R version 4.0.3
getwd()
## [1] "C:/Users/Mufasa/Documents/Repos/MA-climate/code/Rmarkdown"
#setwd("Repos/MA-climate/")
# load data
precip <- readRDS("../../data/processed/decomposed_precip_data_withna.rds")
sst <- readRDS("../../data/processed/decomposed_sst_data_withna.rds")
dim(precip)
## [1] 61200   432
dim(sst)
## [1] 16020   432
# option front or back cut
# for precip front cut
# for sst back cut
cut_matrix <- function(matrix, timelag, side=c("front,back")) {
  if(side == "front") {
    matrix <- matrix[, -c(1:timelag)]
  }
  if(side == "back"){
    matrix <- matrix[, 1:(ncol(matrix)-timelag)]
  }
  return(matrix)
}

# compute correlations
# for different time lags what do we need?
# in: precip path, sst path, time lags
# out: correlation vector
# load precip and sst decomposed
# compute precip means
# remove precip data
# compute correlations with time lags
# correlations out

compute_corr <- function(dec_matrix_sst, dec_matrix_precip, timelag=0) {
  if (timelag!=0) {
    dec_matrix_sst <- cut_matrix(dec_matrix_sst, timelag, "back")
    dec_matrix_precip <- cut_matrix(dec_matrix_precip, timelag, "front")
  }
  assertthat::are_equal(ncol(dec_matrix_precip), ncol(dec_matrix_sst))
  precip_means <- apply(dec_matrix_precip, 2, mean)
  cor_vec <- apply(dec_matrix_sst, 1, function(x) cor(x, precip_means))
  return(cor_vec)
}


#corr0 <- compute_corr(sst, precip,0)
# create matrix/raster
old_sst <-  brick("../../data/interim/sst/ersst_setreftime.nc",
                  varname= "sst")


# plot correlation 
# in: correlation vector, old sst data, timelag for legend
plot_corr <- function(corr_vec, old_sst, timelag,
                      quantiles = c(TRUE,FALSE)) {
  corr_grid <- matrix(corr_vec, nrow = old_sst@nrows, ncol = old_sst@ncols, byrow = TRUE)
  corr_raster <- raster(corr_grid,
                        xmn = old_sst@extent@xmin, xmx = old_sst@extent@xmax,
                        ymn = old_sst@extent@ymin, ymx = old_sst@extent@ymax,
                        crs = old_sst@crs)
  df <- cbind(xyFromCell(corr_raster, 1:ncell(corr_raster)), values(corr_raster))
  df <- base::as.data.frame(df)
  colnames(df) <- c("Longitude","Latitude", "Correlation")
  if(quantiles) {
    q <- quantile(df[,"Correlation"],probs=c(0.025,0.975), na.rm = TRUE)
    df$Correlation[df$Correlation>q["2.5%"] & df$Correlation<q["97.5%"]] <- 0
  }
  plt <- ggplot(data = df, aes(x = Longitude, y = Latitude, fill = Correlation)) +
    annotation_map(map_data("world")) +
    geom_raster(interpolate = TRUE)
  if(!quantiles) {
    plt <- plt #+
      #ggtitle(paste("Correlation of Sea Surface Temperature and Precipitation for Timelag",timelag))
  } 
  if(quantiles) {
    plt <- plt #+
      #ggtitle(paste0("Correlation of Sea Surface Temperature and Precipitation for Timelag"," ",timelag,
                     #",\nValues between quantiles threshold set to 0"))
    
  }
  plt <- plt +
    scale_fill_viridis_c(option="A") +
    theme_bw() +
    coord_quickmap() 
  return(plt)
}

compute_corr_and_plot <- function(
  dec_matrix_sst, dec_matrix_precip, old_sst,
  timelag=0,quantiles = c(TRUE,FALSE)) {
  
  corr_vec <- compute_corr(dec_matrix_sst, dec_matrix_precip, timelag)
  plt <- plot_corr(corr_vec, old_sst, timelag, quantiles)
  return(plt)                                 
}


plot_density <- function(corr_vec, timelag) {
  df <- data.frame(corr_vec)
  plt <- ggplot(df, aes(x=corr_vec)) + geom_density() + ggtitle(paste("Density Plot of Correlations with timelag",timelag))
  plt
}

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
library(jpeg)
library(grid)
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.0.5
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
corr0_vec <- compute_corr(sst, precip,0)
dens0 <- plot_density(corr0_vec, timelag = 0)
ggsave("dens0.jpg", width = 10, height = 10, units = "cm")
## Warning: Removed 5032 rows containing non-finite values (stat_density).
#ggsave("plots/dens0.jpg",dens0)
corr0 <- plot_corr(corr0_vec, old_sst, timelag = 0, quantiles = FALSE)
#ggsave("corr0.jpg", width = 30, height = 30, units = "cm")
ggsave("corr0.jpg", scale = 2)
## Saving 14 x 10 in image
corr0q <- plot_corr(corr0_vec, old_sst, timelag = 0, quantiles = TRUE)
ggsave("corr0q.jpg", scale = 2)
## Saving 14 x 10 in image
# 
# im1 <- rasterGrob(as.raster(readJPEG("dens0.jpg")))
# im2 <- rasterGrob(as.raster(readJPEG("corr0.jpg")))
# grid.arrange(im1,im2, ncol=1)

#ggsave("plots/corr0q.jpg", corr0q)
corr3_vec <- compute_corr(sst, precip, 3)
dens3 <- plot_density(corr3_vec, timelag = 3)
ggsave("dens3.jpg", width = 10, height = 10, units = "cm")
## Warning: Removed 5032 rows containing non-finite values (stat_density).
corr3 <- plot_corr(corr3_vec,old_sst,3,FALSE)
ggsave("corr3.jpg", scale = 2)
## Saving 14 x 10 in image
corr3q <- plot_corr(corr3_vec,old_sst,3,TRUE)
ggsave("corr3q.jpg")
## Saving 7 x 5 in image
# load and save all densities
# dens0
# dens3
corr6_vec <- compute_corr(sst, precip, 6)
dens6 <- plot_density(corr6_vec, timelag = 6)
ggsave("dens6.jpg", width = 10, height = 10, units = "cm")
## Warning: Removed 5032 rows containing non-finite values (stat_density).
corr9_vec <- compute_corr(sst, precip, 9)
dens9 <- plot_density(corr9_vec, timelag = 9)
ggsave("dens9.jpg", width = 10, height = 10, units = "cm")
## Warning: Removed 5032 rows containing non-finite values (stat_density).
corr12_vec <- compute_corr(sst, precip, 12)
dens12 <- plot_density(corr12_vec, timelag = 12)
ggsave("dens12.jpg", width = 10, height = 10, units = "cm")
## Warning: Removed 5032 rows containing non-finite values (stat_density).
captioncaptioncaptioncaptioncaption

caption 

corr3 <- plot_corr(corr3_vec,old_sst,3,FALSE)
ggsave("corr3.jpg")
## Saving 7 x 5 in image
corr3q <- plot_corr(corr3_vec,old_sst,3,TRUE)
ggsave("corr3q.jpg")
## Saving 7 x 5 in image
corr6 <- compute_corr(sst, precip, 6)
plot_corr(corr6,old_sst,6,FALSE)

ggsave("corr6.jpg")
## Saving 7 x 5 in image
plot_corr(corr6,old_sst,6,TRUE)

ggsave("corr6q.jpg")
## Saving 7 x 5 in image
corr9 <- compute_corr(sst, precip, 9)
plot_corr(corr9,old_sst,9,FALSE)

ggsave("corr9.jpg")
## Saving 7 x 5 in image
plot_corr(corr9,old_sst,9,TRUE)

ggsave("corr9q.jpg")
## Saving 7 x 5 in image
corr12 <- compute_corr(sst, precip, 12)
plot_corr(corr12,old_sst,12,FALSE)

ggsave("corr12.jpg")
## Saving 7 x 5 in image
plot_corr(corr12,old_sst,12,TRUE)

ggsave("corr12q.jpg")
## Saving 7 x 5 in image
captioncaptioncaptioncaptioncaption

caption

captioncaptioncaptioncaptioncaption

caption 

library(gtable)
library(grid)
#https://stackoverflow.com/questions/32780320/consistent-figures-size-with-gridextra-in-rmarkdown-knitr-html
g1 <- dens0
g2 <- corr0
g3 <- corr0q
pl <- lapply(list(g1,g2,g3), ggplotGrob)
g123 <- rbind(pl[[1]], pl[[2]], pl[[3]], size="first")
g123$heights <- unit.pmax(pl[[1]][["heights"]], pl[[2]][["heights"]],
                          pl[[3]][["heights"]])
#g34 <- cbind(pl[[3]], pl[[4]], size="first")
#g34$heights <- unit.pmax(pl[[3]][["heights"]], pl[[4]][["heights"]])
#g1234 <- rbind(g12, g34, size="first")
#g123$widths <- unit.pmax(pl[[1]][["widths"]],pl[[2]][["widths"]],pl[[3]][["widths"]])
grid.newpage()
grid.draw(g123)
library(ggpubr)
ggarrange(dens0, corr0, corr0q, ncol=1,nrow=3)
??ggarrange